CONVERGENCE ANALYSIS OF STRANG SPLITTING FOR 
VLASOV TYPE EQUATIONS 



LUKAS EINKEMMER* AND ALEXANDER OSTERMANN* 

Abstract. A rigorous convergence analysis of the Strang splitting algorithm for Vlasov— type 
equations in the setting of abstract evolution equations is provided. It is shown that under suitable 
assumptions the convergence is of second order in the time step h. As an example, it is verified 
that the Vlasov-Poisson equation in 1+1 dimensions fits into the framework of this analysis. Also, 
numerical experiments for the latter case are presented. 
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1. Introduction. The most fundamental theoretical description of a (collision- 
less) plasma comes from the kinetic equation. This so called Vlasov equation is given 
by (see e.g. pQ) 

d t f(t, x,v)+v V/(t, x, v) + F ■ V v f(t, x, v) = 0, 

where x denotes the position and v the velocity. The function / describes a particle- 
probability distribution in the 3 + 3 dimensional phase space. Since a plasma interacts 
with the electromagnetic field in a non-trivial manner, the Vlasov equation needs to 
be coupled to the electromagnetic field. 

Depending on the application, either the full Vlasov-Maxwell equations or a sim- 
plified model is appropriate. Such models include, for example, the Vlasov-Poisson 
and the gyrokinetic equations. 

Due to the high dimensionality of the equations the most common numerical ap- 
proach are so called particle methods. In this class of methods, the phase space is 
left to be continuous and a (large) number of particles with various starting points 
are advanced in time. This is possible due to the structure of the equations, which 
implies that a single particle evolves along a trajectory given by an ordinary differ- 
ential equation. A number of such methods have been developed, most notably the 
particle-in-cell (PIC) method. Such methods have been extensively used for various 
applications (see e.g. [7]). The PIC scheme gives reasonable results in case where the 
tail of the distribution is negligible. If this is not the case the method suffers from 
numerical noise that only decreases as 1/y/n, where n denotes the number of particles 
(see e.g. [12] or [8]). Motivated by the considerations above, a number of schemes 
employing discretization in phase space have been proposed. A comparison of various 
such methods can be found in [5]. 

Using a time splitting scheme for the Vlasov-Poisson equations was first proposed 
by [2] in 1976. In [T^ the method was extended to the Vlasov-Maxwell equations. In 
both cases, second order Strang splitting (see e.g. [13]) is used to advance the solution 
of the Vlasov equation in time. 

Quite a few convergence results are available for semi-Lagrangian methods that 
employ Strang splitting. For example, in [2], [3] and |17j convergence is shown in the 
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case of the 1+1 dimensional Vlasov-Poisson equations. Furthermore, the convergence 
of a special case of the one-dimensional Vlasov-Maxwell equation in the laser-plasma 
interaction context is investigated in [3]. 

In this paper, we will consider a class of Vlasov-type equations as abstract evo- 
lution equation (i.e. without discretization in space). In this context we will derive 
sufficient conditions such that the Strang splitting algorithm is convergent of order 2. 
We will then verify these conditions for the example of the Vlasov-Poisson equations 
in 1+1 dimensions and present some numerical results. 

2. Setting. We will investigate the following (abstract) initial value problem 

(d t f(t) = (A + B)f(t) 

{ /(0) = /o- ( J 

We assume that A is an (unbounded) linear operator and that B can be written in the 
form Bf — B(f)f, where B(f) is an (unbounded) linear operator. We will consider 
this abstract initial value problem on a finite time interval [0, T]. 



Problem (2.1 1 comprises the Vlasov-Poisson and the Vlasov-Maxwell equations 
for A = —v ■ V and appropriately chosen B as special cases. It is also general enough 
to accommodate the gyrokinetic equations (as stated, for example, in |10jV The 
Vlasov-Poisson equations are considered in more detail in section [4] 

2.1. The Strang splitting algorithm. Let fk ~ f(tk) denote the numerical 



approximation to the solution of (2.1| at time t k = kh with step size h. We assume 
that the differential equations d t f = Af and d t g = B k+1 / 2 g, where B k+ i/ 2 is a first- 
order approximation to the operator B (f(t k + -|)), can be solved quite efficiently 
That this is indeed the case for the Vlasov-Poisson equations is investigated in more 
detail in section 03 

The idea of Strang splitting is to advance the numerical solution by the recursion 
fk+i = Skfk, where the (nonlinear) splitting operator S k is given by 

S k =e^'V iJ <=+ 1 /2 e H (2.2) 

The precise conditions on B k+ i/ 2 are given in section [3| below. Resolving this recur- 
sion, we can compute an approximation to the exact solution at time T by 

/„ - ( J] S k J f = S n -i ■ ■ ■ So/o, (2.3) 
\fe=o / 

where n is an integer chosen together with the step size h such that T = nh. 

2.2. Preliminaries. For the convenience of the reader we collect some well 
known results that are used quite extensively in section [3] 

To bound the remainder term Rk(f) of a Taylor expansion 

f(h) = /(0) + hf'(Q) + ... + J^yf^HO) + Rk(f), 
we will use the integral form 

Rk(f) = f f {k \hs){\ - s) k ~ l da, 
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where k > 1. Note that it is implicitly understood that Rk is a function of h as well. 
However, since we will work mostly with a fixed h, it is convenient to drop it in the 
notation of Rk- For convenience we also define 

Mf) = /■ 

For (unbounded) linear operators it is more convenient to work with the ip func- 
tions instead of the remainder term given above. 

Definition 2.1 (ip functions). Suppose that the linear operator E generates a 
Co semigroup. Then we define 

p (hE) = e hE , 

<p k {hE) = f e^ 1 
Jo 



-r)hE 



(fc-1) 



dr, for k > 1. 



(2.4) 



Since we are merely interested in bounds of such functions, we will never di- 
rectly employ the definition given. Instead we will work exclusively with the following 
recurrence relation. 

Lemma 2.2. The if functions satisfy the recurrence relation 

<Pk(hE) = — + hE(p k+1 {hE) 
and in particular (for I G N) 

l-i , k 

e hE = Y t L E k + h l E l <p l (hE). 
^-^ k\ 

k=0 



Proof. The first relation follows from integration by parts applied to (2.4). The 
second by using <po = and applying the first relation repeatedly. □ 

The ip functions are used to expand the exponential of some linear operator. In 
the sense of the previous Lemma, these functions play the same role for an exponential 
of a linear operator as does the remainder term in Taylor's theorem. 

Suppose that the differential equation /' = G(f) has (for a given initial value) 
a unique solution. In this case we denote the solution at time t with initial value 
/(*o) = Jo by the evolution operator Ec(t — to, /o). 

The Grobner-Alekseev formula (also called the nonlinear variation of constants 
formula) will be employed quite heavily. 

Theorem 2.3 (Grobner-Alekseev formula). Suppose that there exists a unique 
f satisfying 

R(f(t)) 




and that f = G(f) has (for a given initial value) a unique solution. Then it holds 
that 

f(t) = E G (t, f ) + ( d 2 E G (t - r, f(r)) R (/(r)) dr. 
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Proof. For linear (and possibly unbounded) G, this formula is proved in [13j by 
the fundamental theorem of calculus. Here, we prove the extension to nonlinear G. 
Let us assume that u(t) is a solution of u'(t) = G (u(t)). By differentiating 

E G (t-T,E G (T-t,u(t)))=u(t) 

with respect to r we get 

-d x E G (t - t, E g (t - t, u(t))) + d 2 E G (t - r, E g (t - t, u(t))) G {E g (t - t, u(t))) = 0, 
which can be rewritten as 

-diE G (t - t, u(t)) + d 2 E G (t - t, u(t)) G (u(t)) - 0. 
The initial value of u is now chosen such that u(t) = f(r) which implies 

-d t E G (t - t, /(r)) + 8 2 E G (t - r, /(r)) G (/(r)) = 0. 

Altogether we have for (p{r) = E G (t — r, /(r)) (by the fundamental theorem of cal- 
culus) 

f(t)-E G (tJ )= [ ip'(T)dT 



(-diE G (t - r, f(r)) + d 2 E G (t - r, /(r))/'(r)) 
d 2 E G (t-r,f(T))R(f(T)) dr, 



as desired. □ 

Since in some expansions anticommutator relations appear quite naturally, we 
will employ the notation 

{Ei,E 2 } = E X E 2 + E 2 Ei, 

for linear operators E\ and E 2 (on a suitable domain). 

In what follows C will denote a generic constant that may have different values 
at different occurrences. 

3. Convergence analysis in the abstract setting. The problem of splitting 
an evolution equation into two parts, governed by linear and possibly unbounded 
operators, has already been investigated in some detail. In it is shown that 
splitting methods with a given classical order retain this order in the stiff case (under 
suitable regularity assumptions). 

An alternative analysis for Strang splitting in the linear case is given in |14) . The 
approach presented there is more involved, however, it demands less regularity on the 
solution. The purpose of this section is to extend this analysis to the abstract initial 
value problem given by (2.1 1. 

3.1. Convergence. Let us begin by defining a suitable notion of consistency for 
our splitting operator. 

Definition 3.1 (Consistency of order p). The Strang splitting algorithm (2.2 1 
is consistent of order p if 

\\f(t k + h)-S k f(t k )\\x <ChP+\ 
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holds uniformly in < t/- = kh < T. 

Note that for algorithm (2.2 1, the order of consistency is not necessarily p = 2. 
The actual order depends on the properties of the involved operators, and order 
reduction can show up even in the linear case, see [14] . 

Theorem 3.2 (Convergence). Suppose that the Strang splitting (2.2) is consistent 
of order p and non-expansive, i.e. 



\Sk\\ 



< 1. 



(3.1) 



Then it is convergent of order p, i.e. ||/„ — f{t n )\\x < Ch p . 

Proof. The proof is quite standard. We rewrite the global error by means of the 
following telescopic identity 



H Sfc fo - f(nh) = ]T J] S n ( S kf{hk) - f(hk + h)) . 



Then from the consistency bound and by assumption ( |3.1[ ) (which implies stability) 
we follow 



fn-l 



n s * /o - f( nh ) 



\k=0 



which is the desired bound. □ 



n-l 



<J2\\(S k f(hk)-f(hk + h))\\ x < Ch 2 
x fe=o 



3.2. Consistency. It is the purpose of this section to formulate assumptions un- 



der which the consistency bound holds for the abstract initial value problem (2.1 ). To 



make the derivations less tedious we will adhere to the notation laid out in Remark l3.31 
Remark 3.3. In this section we will denote the solution of (2.1 1 at a fixed time 
tk by /o- The notation f(r) is then understood to mean f(t k + r). The function 
fo is a (possible) initial value for a single time step (i.e., a single application of the 
splitting operator Sk). It is not, in general, the initial value of the solution to the 
abstract initial value problem as in the previous sections. If we assert that a bound 
holds uniformly in tk, it is implied that it holds for all fo in the sense defined here 
(remember that tk £ [0, T]). Since tk is fixed we will use the notation -B1/2 and S 
instead of Bk+1/2 and Sk respectively. 

Let us start with expanding the exact solution by using the Grobner-Alekseev 
formula (this has been proposed in the context of the nonlinear Schrodinger equation 
in [15] ) . It is of emphasize that we consider the linear operator A as a perturbation of 
the differential equation given by the non-linear operator B. This choice is essential 
for the treatment given here, since it allows us to apply the expansion sequentially 
without any additional difficulties. 

Lemma 3.4 (Expansion of the exact solution). 
the formal expansion 



The exact solution of (2.1) has 



f(h) = E B (h,f )+ / d 2 E B (h-r,f(T))AE B (rJo)dT 
Jo 

rh i-t 



+ 



[ [ d 2 E B (h-T,f(T))Ad 2 E B (T-aJ(a))AE B (a,f )dadT 

Jo Jo 

I I 1 I \ fl d ^ E B(n - 7*+i, f{r k+ x))A ] /(r 3 ) dT 3 dr 2 dn, 

Jo Jo Jo V fe=0 ) 
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where tq = h. 

Proof. Apply the Grobner-Alekseev formula three times to equation (2.1 1. □ 
Next we expand the splitting operator S in a form that is suitable for comparison 

with the exact solution. 

Lemma 3.5 (Expansion of the splitting operator). The splitting operator S has 

the formal expansion 

Sfo = e hB ^ h + \{A, e hB ^ } f + ^ { A, { A, e hB ^ }}f + R 3 fo, 



where 



16 Jo 



(1-s) 2 {a,{a,{ 



A,e^ A e hBl ' 2 e^ L 



*.}}} 



As. 



Proof. Let us define the function g(s) = Q2 sA e hB ^/2 e l sA . The first three deriva- 
tives of g are given by 

g'(s) = l{A,g(s)}, 
g"( S ) = ±{A,{Ag(s)}}, 
9^(s) = U A dA{Ag(s)}}}. 



From the observation that S = g(h) and by Taylor's theorem we obtain the result. □ 
In the remainder of this section, we use the notation 

R i (B)g = h- i (B-B 1/2 )g, ie{0,l,2}. 

Theorem 3.6 (Consistency). Suppose that the estimates 

\\A i ^ h - s ) B ^R % _ i {B)E B {s,f )\\x<C, ie {0,1,2} (3.2) 
||(%i2i(B))/o|U < C, (3.3) 
\\A^B)+^ Vl+5M {hB 1/2 )A 1+ ^f \\ x <C, ie {0, 1,2} (3.4) 

\\A^R 1+Sw (d 2 E B (;f Q ))A 1+ ^f \\ x <C\ ie {0,1,2} (3.5) 

hold uniformly in t and in s € [0, ft]. In addition, suppose that for all kj G N, with 
Siii — ^ — h the estimates 



l{Dfd 2 E B (s J ,f(a J ))A d^E B (s i+l ,f ) 



v =1 



<C, ie{l,2} (3.6) 



A" 



Y[d 2 E B (s k -a k J(a k ))A) f(s) 



vfc=l 



<c, 



x 



{A,{A,{A,et A e hB ^ei A }}}f \\ x <a 



(3.7) 
(3.8) 



hold uniformly in t as well as in s € [0, ft], Sj € [0, ft], and o~j € [0, ft], where D^ 3 is a 
differential operator of order kj in the variables Sj and o~j . Then the Strang splitting 
(2.2) is consistent of order 2. 
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Proof. We have to compare terms of order 0, 1, and 2 in Lemma [3.4| and Lemma [3.5| 
and show that the remaining terms of order 3 can be bounded as well. 

Terms of order 0. We can rewrite E b (t, fo) as u(t), where u(t) satisfies the 
differential equation 

v! = B 1/2 u + (B - B 1 / 2 )u, 
with initial value fo- Employing the variation of constants formula we get 



u(h) 



7 



o+ / e^ B ^(B-B 1/2 )E B (T,f )dr. 
Jo 



Assumption (3.2) with i = gives the desired order. 
Terms of order 1 . For 

? (r)=e( fc -^'/Me^./o, k(r)=d 2 E B (h-rJ(r))AE B (r,f ) 

we get (by use of the trapezoidal rule) 

h 



Hr) dr=^ b(0) - fc(0) + g(h) k(h)]~ 



6{l-0)k"{6h) d6. 



First, let us compare g(h) and k(h) 

g(h) - k(h) = A [e' iS V2 /o _ E B (h, f )] , 

which is the same as before, except that we have an additional A to the left of the 
expression. We are left with assumption (3.2) with i = 1. Second, we have to compare 
g(0) and fc(0) 

g(Q) - k(0) = d 2 [e hB ^ _ E B (h, f )] Af , 

where we can pull out <9 2 as it does not effect the linear operator e hBl / 2 . To handle 
the derivative with respect to the initial value, we expand both terms (using the if 
functions for the exponential) and get 

, 9 (0) - fc(0) = h 2 Ud2Ri(B))f + B 2 1/2 f 2 (hB 1/2 ) - R 2 (d 2 E B (-, /„))] Af . 



The first term is bounded by assumption (3.3). The second term is bounded by 
assumption (3.4) with i = and the third term by assumption (3.5) with i = 0. 

Finally, we have to estimate the remainder term of the quadrature rule which is 
bounded by assumption (3.6) with i = 1. 

Terms of order 2. For the functions g(r, a) = e {h ^ T ^ Bl ' 2 Ae {T ~^ Bl ' 2 Ae aBl ' 2 f and 
k(r, a) = d 2 E B (h — r, f(T))Ad 2 E B (r — a, f(a))AE B (a, fo) we employ a quadrature 
rule 

h r r 



y [g(0, 0) + 2g(h, 0) + g(h, h)] - J J fc(r, a) dadr 



\g(0, 0) + 2g(h, 0) + g(h, h) - fe(0, 0) - 2k(h, 0) - k(h, h)] + d, 



8 



L. EINKEMMER AND A. OSTERMANN 



where d is the remainder term. From this, it can be shown that we have to bound 
the following three terms 

A 2 [e hB ^fo-E B (h,f Q )] 
[e hB v* - d 2 E B {h, /„)] A 2 f = h [B 1/2 ^(hB 1/2 ) - Rx{d 2 E B {-, / ))] A 2 f 
2A [e hB ^ - d 2 E B (hJ )] Af = 2hA [B 1/2Vl (hB 1/2 ) - R 1 (d 2 E B (-, f ))] Af . 



The first term we can bound by using assumption (3.2) with i — 2. In addition 



we can bound the second and third term using assumption (3.4) with i = 1 and i = 2 



and assumption |3.5| with i = 1 and i = 2 respectively. Finally, the remainder term 



depends on the partial derivatives of k(r,a) and can be bound by (3.6) with i = 2 



Terms of order 3. In order to bound the remainder terms in the expansion of 



the exact solution as well as the approximate solution, we need assumption (3.7) and 



(3.8) respectively. □ 



4. Convergence analysis for the Vlasov Poisson equations. We will con- 
sider the Vlasov-Poisson equations in 1+1 dimensions, i.e. 

d t f(t, x, v) = -vd x f(t, x, v) - E(f(t, ■, -),x)d v f(t, x, v) 

d x E(f(t,-,-),x)=Jj(t,x,v)dv-l 

f(0,x,v) = fo(x,v), 

with periodic boundary conditions in space. The domain of interest has length L. 
Thus, for all x G R 

f(t,x,v) = f(t,x + L,v). 

As will be apparent in the next section it is unnecessary to impose boundary conditions 
in the velocity direction. This is due to the fact that for a function f with compact 
support in the velocity direction the solution will continue to have compact support 



for all finite time intervals [0, T] (see Theorem 4.1) 



4.1. Definitions and notation. The purpose of this section is to introduce the 
notations and mathematical spaces necessary for giving existence, uniqueness, and 
regularity results as well as to conduct the estimates necessary for showing consistency 
and stability. 

For estimation we will use the Banach space L 1 ([0, L] x R) exclusively. This is 
reasonable as the solution / represents a probability density function. As such the L 1 
norm is conserved for the exact (as well as the approximate) solution. Nevertheless, 
all the estimations could be done as well, for example, in L°°([0, L] x R). 

However, we need some regularity of the solution. This can be seen from the 
assumptions of Theorem |3.6[ where we have to apply a number of differential oper- 
ators to the solution f(t). Thus, we introduce the following spaces of continuously 
differentiable functions 

C™r,c = {/ G C m (R 2 , R) ; Va, v. (f(x + L,v) = f(x, v)) A (supp/fo ■) compact)} , 
C£r = {/ e C" l (R,R) ; Vx: f(x + L) = f(x)} . 

Together with the norm of uniform convergence of all derivatives up to order m, the 
spaces C™ r c and C pcr are turned into Banach spaces. 
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We also have to consider spaces that involve time. To that end let us define 

C™(0,T;C m ) = \feC m ([0,T\,C°); (f(t)eC m )A( sup ||/(t)|| Cm < oo) 1 , 
{ te[o,T] J 

where C m is taken as either C™ r c or C™ r . It should be noted that if it can be shown 
that the solution / of the Vlasov-Poisson equations lies in the space C m (0,T;C m ), 
we can bound all derivatives (in space) up to order m uniformly in t 6 [0, T]. 

4.2. Existence, uniqueness, and regularity. In this section we recall the 
existence, uniqueness, and regularity results of the Vlasov-Poisson equations in 1+1 
dimensions. The theorem is stated with a slightly different notation in [3] and [2]. 

Theorem 4.1. Assume that f € C™ r c is non-negative, then f £ C m (0, T; C™ r c ) 
and E(f(t,-,-),x) as a function of(t,x) lies in C m (0, T; Cp™ r ) . /n addition, we can find 
a number Q(T) > smc/i that for all t € [0, T] and :r G R ZioWs that supp/(t, x, •) C 
[-Q(T),Q(T)}. 

Proof. A proof can be found in [§]. □ 

We also need a regularity result for the electric field that does not directly result 
from a solution of the Vlasov-Poisson equations, but from some generic function / 
(e.g., computed from an application of a splitting operator to /o). 
Corollary 4.2. For f e C™ t . c it holds that E(f, •) € C^ r . 



Proof. The result follows from the proof of Theorem 4.1 In addition, in the 



1+1 dimensional case it can also be followed from the exact representation of the 



electromagnetic field that is given in (5.2) below. □ 



It should also be noted that due to the proof of Theorem |4.1| the regularity results 
given can be extended to the differential equations generated by B and B 1 / 2 - Thus, 
remains ,/o) or e tBl ' 2 f is substituted for f(t). 



Theorem 



4.1 



4.3. Consistency. The most challenging task in proving the assumptions of 



Theorem 3.6 is to control the derivative of Eb with respect to the initial value. The 
following lemma accomplishes that. 
Lemma 4.3. The function 

nra v pt . /?min(m— 1,-tj 

*^per,c ^pcr.c * ^pcr,c 

(uo,g) >-> d 2 E B (t,u )g, 

is well defined. 

Proof. We solve u'(t) = Bu(t) and it(0) = uq. Motivated by the method of 
characteristics we can write 

Vl (t) = E(u(t,;-),x) 

Vu o (0)=v 
u(t,x,v) = u (x,V Uo (t)(x,v)). 

We now compute the Gateaux derivative with respect to the direction g and get 

d h E B (t,u + hg)(x,v)\h=o = {d 2 u ) (x, V Uo (t)(x,v)) V g (t)(x,v) + g(x, V u „(t)(x, v)), 

since V is linear with respect to the initial value. From this representation the result 
follows. □ 

In addition, it is usefuel to have an explicit representation of the Taylor series 
remainder terms of B and Bi/ 2 - However, up until now we have only specified that 
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Bx/2 is a first order approximation to B. To perform the analysis that follows we will 
approximate B 1 / 2 by an evaluation of B(-). The details are given in the following 
definition. 

Definition 4.4 (Structure of B 1/2 ). Let us define 

B 1/2 .f = B(f h/2 )f, 

where fh/ 2 is an approximation to f(h/2) of order 1. 

It remains, however, to verify that this ansatz produces an approximation of 
order 1, i.e. it is consistent with the assumptions we made above. 

LEMMA 4.5. B 1 / 2 is an approximation of order 1 to B. 

Proof. Since f(h/2) = /(0) + |/'(0) + O (h 2 ) the assumption implies that f h/2 = 
/(0) + f /'(0) + O (h 2 ). Expanding B we get 

B (/ (J)) = S(/(0)) + ^'(/(0))/'(0) + O (h 2 ) . 
Comparing this with the expansion of B 1 / 2 

B(f h/2 ) = B(/(0)) + £i?'(/(0))/'(0) + O (h 2 ) 
our result follows. □ 

The following two lemmas present time derivatives up to order 2 of Bf, Bf h / 2 
and E B (t, f ) which follow from a simple calculation. Let us start with the derivatives 
of the operator B and B 1 / 2 applied to some f(t) = f(t, •, •). 

Lemma 4.6. For f sufficiently often continuously differentiable, we have 

d t Bf(t, x, v) = E (f(t), x) d v f(t, x, v) + E(f(t),x)d v f(t, x, v) 

d 2 t Bf(t, x, v) = E (f"(t),x) d v f(t, x, v) 

+ 2E (f'(t),x) d v f'{t, x, v) + E(f(t),x)d v f"(t, x, v) 
d t B 1/2 f(t,x,v) = E(f h/2 ,x)d v f(t,x,v) 
d 2 B 1/2 f{t, x, v) - E(f h/2 ,x)d v f"(t, x, v) 



Proof. From the relation Bf(t, x, v) = E(f(t), x)d v f(t, x, v)) and B 1 / 2 f(t, x, v) — 
E(fh/ 2 ,x)d v f(t,x,v) the result follows by the product rule. □ 

Also we have to compute some derivatives of the evolution operator Es{t, fo) 
with respect to time. 

Lemma 4.7. For f sufficiently often continuously differentiable, we have 

d t E B (t,f )=BE B (t,f ) 

d 2 E B {t, f ) = E(E B (t, f ),x)d v (BE B (t, f )) + E(BE B (t, f ),x)8 v E B (t, fo) 
dt(d 2 E B (tJ )) - E(E B (t,fo),x)d v (d 2 E B (t,f ))+E(d 2 E B (t,f ),x)d v E B (t,f ) 
d 2 t {d 2 E B {t,fo)) = E(BE B (tJo),x)d v (d 2 E B (tJ )) 

+ E(E B (t, fo), x)d v (d t (d 2 E B (t, f ))) 

+ E(d t (d 2 E B (t, f )),x)d v E B (t, f ) 

+ E(d 2 E B (t,f ),x)d v BE B (t,fo) 
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Proof. From the relation Bf(t, x, v) — E(f(t), x)d v f(t, x, v)) the result follows by 
a simple calculation. □ 

It is also necessary to investigate the behavior of the ip functions introduced in 
Definition |2~T1 

Lemma 4.8. For the Vlasov-Poisson equations the functions E min ^ h1 ^ <pi(hE) 
with E € {A, -B1/2} are maps from C™ r c to C™ r c for all i € N. 
Proof. For i = we have 



-hvd a 



fo(x,v) = f (x - vh,v), 



and 



-hE(f h/ M foM = /q ( XjV _ E (f h/2 ,x)h) 



This clearly doesn't change the differentiability of /q. For the ip functions we proceed 



by induction. Since by Definition 2.1 we have 



hEip k+1 (hE) = <p k (hE) 



I. 



from which the desired result follows. □ 

Now we are able to show that all the assumptions of Theorem [3^6] are fulfilled and 
that we thus have consistency of order 2. This is the content of theorem |4.9| 

Theorem 4.9. Suppose fo € Cp Crc , fo is non-negative and f h / 2 is an approxi- 
mation to f(h/2) of order 1. Then Strang splitting for the Vlasov-Poisson equations 
is consistent of order 2. 

Proof. The proof proceeds by noting that the solution has compact support (for 
a finite time interval), i.e. we can estimate v by some constant Q. On the other 
hand it is clear that for f e C™+* we get Af £ C™ r c and B 1/2 fo € C™ r c . The 

terms 
the tp 



4.6 



4.8 



same is true for B as can be seen by Corollary |4.2| Noting that by Lemma 
of the form Ri(B) are mappings from C™+* to C™ r c and that by Lemma 
functions are mappings from Cp" r c — » C™ T c we can conclude that after applying all 
operators in assumptions (3.2), (3.3), (3.4), and (3.5) we get a continuous function. 



By the regularity results we can bound these functions uniformly in time. The same 



argument also shows the validity of the bound in assumption (3.8) 



Finally, with the help of Lemma [4~3| an d [4~7| together with the above observations 
we can bound assumptions (3.6) and (13. 71). □ 



4.4. Convergence. We are now in the position to prove second order conver- 
gence of Strang splitting for the Vlasov-Poisson equations in L 1 . The same result 
holds literally in L°° (or any other LP space). 

Theorem 4.10. Suppose fo € Cp Crc , fo is non-negative and fhji is an approxi- 
mation to f(h/2) of order 1. Then Strang splitting for the Vlasov-Poisson equations 
is convergent of order 2. 

Proof. By Theorem |4.9| and |3.2| it is sufficient to show that 



gl j4 g' l - B fc + l/2 e f A j 



< 11/11 



It is evident from the proof of Lemma |4.8| that the above operators can be represented 
as a translation only. Since a translation does not change the L 1 norm our result 
follows. □ 
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5. Numerical experiments. In this section we present some numerical exper- 
iments. Even if we neglect space discretization for the moment, we have to settle on 
some fh/2 that is an approximation of f(h/2) of order 1. This can be achieved by 
Taylor series expansion, interpolation of previously computed values, or by making an 
additional Lie-Trotter time step of length h/2. Since we are interested in time integra- 
tion only, we choose the latter. This method is trivial to implement (once the Strang 
splitting scheme is implemented) and doesn't suffer from the numeric differentiation 
problems of a Taylor expansion. Thus, one possible choice is to use 

A/2 = e* fl <'°>e^/o (5.1) 
in our simulations. However, since the semigroup generated by B(/q) can be repre- 



sented as a translation in velocity (see the proof of Lemma 4.8 1 and the electric field 



depends only on the average of the density function with respect to velocity (i.e. it 



depends only on the charge density), it is possible to drop the first factor in (5.1| 



without affecting the resulting electric field. Since the computation of the second 



operator in (5.1) is the first step in the Strang splitting algorithm, this leads to a 
computationally efficient scheme. This scheme is also employed in [16] . for example. 
However, no argument why second order accuracy is retained is given there. 
To compute the electric field we will use the following formula (see e.g. [3]) 

E(f(t, -,-),x)= [ K(x, y)( [ f(t, y, v)dv - l) dy, 




1 n , ( 5 ' 2 ) 
1 < x < y, 

y < x < L. 

For space discretization we will employ the discontinuous Galerkin method described 
in [IB]. The approximation is of order 2 with 80 cells in both the space and velocity 
direction. In [16 the coefficients up to order 2 are given. However, it is no difficult 
to use a computer program to compute the coefficients to arbitrary order. 

5.1. Landau damping. The Vlasov-Poisson equations in 1+1 dimensions to- 
gether with the initial value 

fo{x,v) = — Le- l,2/2 (l + acos(0.5a;)), 
V 2vr 

is called Landau damping. For a = 0.01 the problem is called linear or weak Landau 
damping and for a = 0.5 it is referred to as strong or non-linear Landau damping. As 
can be seen, for example, in [H1IS] and [18] Landau damping is a popular test problem 
for Vlasov codes. 

For comparison we display the error of the Strang splitting algorithm together 
with the error for Lie-Trotter splitting (a method that is expected to be of order 1). 
Since we are mainly interested in the time integration error and there is no analytical 
solution of the problem available, we compare the error for different step sizes with a 
reference solution computed with h = 3.9 • 10~ 3 . The error is computed in the discrete 



L 1 norm at time t = 1. The results are given in Figure 5.1 



6. Conclusion. In this paper sufficient conditions are given that guarantee con- 
vergence of order 2 for the Strang splitting algorithm in the case of Vlasov-type 
equations. It is then shown that the Vlasov-Poisson equations in 1+1 dimensions 
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Fig. 5.1. Error of the particle density function /(l,-,-) /or Strang and Lie-Trotter splitting 
respectively, where a = 0.01 (top) and a = 0.5 (bottom). 



is an example of a Vlasov-type equation, i.e. that it fits into the framework of the 
analysis conducted. For simulation on a computer, however, an approximation has 
to be made (i.e. some sort of space discretization has to be introduced). This ap- 
proximation is not included in the analysis done here. Nevertheless, the numerical 
experiments suggest that second order convergence is retained in the fully discretized 
case as well. 
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